Skip to contents

happign for foresters

Source: vignettes/web_only/happign_for_foresters.Rmd
happign_for_foresters.Rmd

First choose a zone of interest

For the example we will work with the Gâvre forest. The starting point come from google map : you just have to select a random point inside the forest by using right click on map.

point = st_sfc(st_point(c(-1.7977020162748798, 47.52634370369097)), crs = 4326)
tm_shape(point)+
   tm_dots(size = 0.1)+
   tm_basemap("OpenStreetMap")

Extract commune borders

Now that we have point, getting shape of commune is just a matter of finding the good layer name from IGN. The one we are after is called “ADMINEXPRESS-COG-CARTO.LATEST:commune” from the “Administrative” category. Then we use the get_wfs() function to download the shape.

As you can see below, there are three communes close to our point as you can see below (tip : click on the interactive map below to get population). For other example we will only select one of them.

apikey <- get_apikeys()[1]  # "administratif"
layer_name <- "ADMINEXPRESS-COG-CARTO.LATEST:commune"

# The layer_name I use come from `all_layers = get_layers_metadata(apikey, "wfs") that return all of them layer_name

borders <- get_wfs(shape = point,
                   apikey = apikey,
                   layer_name = layer_name)
#> Request 1/1 downloading...

tm_shape(borders)+
   tm_polygons(lwd = 2,
               col = "nom",
               id = "nom",
               popup.vars = "population")+
   tm_basemap("OpenStreetMap")

Download cadastral parcel

The cadastral parcel is the basis of any forest management, it is essential. With happign, it is easy to download it :

# Select only commune to reduce downloading time for hight resolution raster
borders <- borders %>% 
   slice(1)

apikey <- get_apikeys()[14] # "parcellaire"
layer_name <- get_layers_metadata(apikey, "wfs")
name_parcellaire_layer <- layer_name[15,2] # "CADASTRALPARCELS.PARCELLAIRE_EXPRESS:parcelle"

parcellaire <- get_wfs(shape = borders,
                       apikey = apikey,
                       layer_name = name_parcellaire_layer) %>% 
   st_transform(st_crs(borders)) %>% 
   st_intersection(borders)
#> Request 1/10 downloading...
#> Request 2/10 downloading...
#> Request 3/10 downloading...
#> Request 4/10 downloading...
#> Request 5/10 downloading...
#> Request 6/10 downloading...
#> Request 7/10 downloading...
#> Request 8/10 downloading...
#> Request 9/10 downloading...
#> Request 10/10 downloading...
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

tm_shape(borders)+
   tm_borders(col = "red", lwd = 2)+
tm_shape(parcellaire)+
   tm_polygons(alpha = 0)+
   tm_basemap("OpenStreetMap")

Downloading BD Forêt

The first interesting layer for forester is the “BD Forêt” which is all vegetation type assigned to each area greater than or equal to 0.5 ha (5,000 m²). There is two layer for forest : the old one BD Forêt V1 and the new one BD Forêt V2.

apikey <- get_apikeys()[9]  #"environement"
layer_name <- get_layers_metadata(apikey, "wfs")
name_BDV1 <- layer_name[1,2] #LANDCOVER.FORESTINVENTORY.V1:resu_bdv1_shape
name_BDV2 <- layer_name[2,2] #LANDCOVER.FORESTINVENTORY.V2:bdforetv2

BDV1 <- get_wfs(borders, apikey, name_BDV1) %>% 
   st_transform(st_crs(borders)) %>% 
   st_intersection(borders) # Layer is download for a bbox around the shape. An intersection is done to get forest only inside our zone of interest
#> Request 1/1 downloading...

BDV2 <- get_wfs(borders, apikey, name_BDV2) %>%
   st_transform(st_crs(borders)) %>% 
   st_intersection(borders)
#> Request 1/1 downloading...

tm_shape(BDV1) +
   tm_polygons(col = "libelle",
               popup.vars = names(BDV1)[1:20],
               legend.show = FALSE) +
tm_shape(BDV2) +
   tm_polygons(col = "tfv",
               alpha = 0.5,
               popup.vars = names(BDV2)[1:16],
               legend.show = FALSE) +
tm_shape(borders) +
   tm_borders(lwd = 2) +
   tm_basemap("OpenStreetMap")

With the map below, you can choose which BD Forêt you want to display. The BDV2 is obviously more precise. Many differences can be calculated, one is the area : BDV2 as 285719.091648482 more hectare described.

More calculations can be done as you can see below :

forest_type_BDV2 = BDV2 %>%
  mutate(area = as.numeric(st_area(geometry))) %>%
  st_drop_geometry() %>%
  group_by(essence) %>%
  summarise(sum_area = sum(area)/10000) %>%
  arrange(sum_area) %>%
  mutate(essence = as.factor(essence))

ggplot()+
  geom_col(data = forest_type_BDV2,
           aes(x = rev(reorder(essence, sum_area)),
               y = sum_area,
               fill = as.factor(essence)))+
  theme_bw()+
  labs(title = "Surface couverte par essences [ha]",
       y = "Surface [ha]",
       fill = "Essence :")+
  theme(axis.text.x = element_blank())

Detect protected area

One information you really want when you work at forest management is if your zone of interest is inside protected area. The example code below is design to test automatically test every layer starting with “PROTECTED” so you’re can be sure that you have all of them.

So, for the Brocéliande forest we find : - 4 forest rescue point - 1 SIC ("Sites d’Importance Communautaire) - 15 Znieff 1 (blue on map) - 1 Znieff 2 (grey on map)

Again, you can click on map, point and shape for more informations.

apikey <- "environnement"
protected_area_names <- get_layers_metadata(apikey, "wfs") %>% 
   filter(grepl("^PROTECTED", name)) %>% 
   pull(name)

all_protected_area <- map(.x = protected_area_names,
                         .f = ~ try(get_wfs(borders, apikey, .x), silent = TRUE)) %>% 
   set_names(protected_area_names) %>% 
   discard(~ is.null(dim(.)))
#> Request 1/1 downloading...
#> Request 1/1 downloading...
#> Request 1/1 downloading...
#> Request 1/1 downloading...
#> Request 1/1 downloading...

# Plot the result
tm_shape(all_protected_area[[1]])+
   tm_dots(group = "Point rencontre des secours en forêts", col = "red")+
tm_shape(all_protected_area[[2]])+
   tm_polygons(group = "SIC", alpha = 0.8)+
tm_shape(all_protected_area[[4]])+
   tm_polygons(group = "Znieff 2", alpha = 0.8)+
tm_shape(all_protected_area[[3]], alpha = 0.8)+
   tm_polygons(group = "Znieff 1", alpha = 0.8, col = "blue")+
tm_shape(borders) +
   tm_borders(lwd = 2) +
   tm_basemap("OpenStreetMap")

MNS, MNT and MNH…

It’s always good to know a more about the terrain topologie. The ign offers a MNT and a MNS for download. As a reminder, the MNT corresponds to the surface of the ground and the MNS to the real surface (in our case, the trees). It is thus easy to find the height of the trees by subtracting the DTM from the MNS.

layers_name <- get_layers_metadata("altimetrie", "wms")
MNT_name <- layers_name[2,2] # ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES
MNS_name <- layers_name[3,2] # "ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES.MNS"

MNT <- get_wms_raster(borders, "altimetrie", resolution = 5, MNT_name, filename = "MNT")
#> 1/2 downloaded
#> 2/2 downloaded
file.remove("MNT.tif") # I don't want to keep MNT on disk because this is vignette
#> [1] TRUE
MNS <- get_wms_raster(borders, "altimetrie",resolution = 5, MNS_name, filename = "MNS")
#> 1/2 downloaded
#> 2/2 downloaded
file.remove("MNS.tif") # I don't want to keep MNS on disk because this is vignette
#> [1] TRUE

MNH <- MNS - MNT
MNH[MNH < 0] <- NA  # Remove negative value 
MNH[MNH > 40] <- NA # Remove height more than 40m

tm_shape(MNH)+
   tm_raster()+
tm_shape(borders)+
   tm_borders(lwd = 2)+
   tm_basemap("OpenStreetMap")
#> stars object downsampled to 980 by 1020 cells. See tm_shape manual (argument raster.downsample)

NDVI

The code below present the calculation of the NDVI. All informations and palette come from this website The value range of an NDVI is -1 to 1. It is (Near Infrared - Red) / (Near Infrared + Red) :

  • Water has a low reflectance in red, but almost no NIR (near infrared) reflectance. So the difference will be small and negative, and the sum will be small, and NDVI large and negative.
  • Plants have a low reflectance in red, and a strong NIR reflectance. So the difference will be large and positive, and the sum will be just about the same as the difference, so NDVI will be large and positive.

Categories are somewhat arbitrary, and you can find various rules of thumb, such as:

  • Negative values of NDVI (values approaching -1) correspond to water. Values close to zero (-0.1 to 0.1) generally correspond to barren areas of rock, sand, or snow. Low, positive values represent shrub and grassland (approximately 0.2 to 0.4), while high values indicate temperate and tropical rainforests (values approaching 1).
  • Very low values of NDVI (0.1 and below) correspond to water, barren areas of rock, sand, or snow. Moderate values represent shrub and grassland (0.2 to 0.3), while high values indicate temperate and tropical rainforests (0.6 to 0.8).
# To show the possibility of 20cm resolution for IRC, let's take only the biggest parcels
biggest_parcels <- parcellaire %>% 
   mutate(area = st_area(.)) %>% 
   slice_max(area)

IRC <- get_wms_raster(shape = biggest_parcels,
                      apikey = "ortho",
                      resolution = 0.2,
                      layer_name = "ORTHOIMAGERY.ORTHOPHOTOS.IRC",
                      filename = "IRC")
#> 1/4 downloaded
#> 2/4 downloaded
#> 3/4 downloaded
#> 4/4 downloaded
file.remove("IRC.tif") # I don't want to keep IRC on disk because this is vignette
#> [1] TRUE

infrared = IRC %>% slice(band, 1)
red = IRC %>% slice(band, 2)

NDVI = (infrared-red)/(infrared+red)
   
breaks_NDVI = c(-1,-0.2,-0.1,0,0.025 ,0.05,0.075,0.1,0.125,0.15,0.175,0.2 ,0.25 ,0.3 ,0.35,0.4,0.45,0.5,0.55,0.6,1)
   
palette_NDVI =  c("#BFBFBF","#DBDBDB","#FFFFE0","#FFFACC","#EDE8B5","#DED99C","#CCC782","#BDB86B","#B0C261","#A3CC59","#91BF52","#80B347","#70A340","#619636","#4F8A2E","#407D24","#306E1C","#216112","#0F540A","#004500")

tm_shape(borders)+
   tm_borders(lwd = 2, col = "red")+
tm_shape(NDVI)+
   tm_raster(stretch.palette = F,
             breaks = breaks_NDVI,
             palette = palette_NDVI,
             colorNA = "red")+
tm_shape(biggest_parcels)+
   tm_borders(lwd = 2, col = "blue")+
   tm_basemap("OpenStreetMap")
#> stars object downsampled to 1031 by 970 cells. See tm_shape manual (argument raster.downsample)

The gloss index

The gloss index represents the average of the image glosses. This index is therefore sensitive to the brightness of the soil, related to its moisture and the presence of salts on the surface. It characterizes especially the albedo (solar radiation that is reflected back to the atmosphere). The gloss index allows us to estimate whether the observed surface feature is light or dark.


gloss_index <- sqrt(red^2 + infrared^2)

tm_shape(borders)+
   tm_borders(lwd = 2, col = "red")+
tm_shape(gloss_index)+
   tm_raster()+
tm_shape(biggest_parcels)+
   tm_borders(lwd = 2, col = "blue")+
   tm_basemap("OpenStreetMap")
#> stars object downsampled to 1031 by 970 cells. See tm_shape manual (argument raster.downsample)

Last but not least… BD Topo

BD topo from IGN covers in a coherent way all the geographical and administrative entities of the national territory. So you can find in it :

  • Administrative (boundaries and administrative units);
  • Addresses (mailing addresses) ;
  • Building (constructions) ;
  • Hydrography (water-related features) ;
  • Named places (place or locality with a toponym describing a natural space or inhabited place);
  • Land use (vegetation, foreshore, hedge);
  • Services and activities (utilities, energy storage and transportation, industrial sites);
  • Transportation (road, rail or air infrastructure, routes);
  • Regulated areas (most of the areas are subject to specific regulations).

For the example below I choose to eownload all water-related data :

cour_eau <- get_wfs(borders, "topographie", "BDTOPO_V3:cours_d_eau") %>% 
   st_intersection(st_buffer(borders,100))
#> Request 1/1 downloading...
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
detail_hydro <- get_wfs(borders, "topographie", "BDTOPO_V3:detail_hydrographique") %>% 
   st_intersection(st_buffer(borders,100))
#> Request 1/1 downloading...
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

# water detected by satellite
surf_hydro <- get_wfs(borders, "topographie", "BDTOPO_V3:surface_hydrographique") %>% 
   st_intersection(st_buffer(borders,100))
#> Request 1/1 downloading...
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

tm_shape(cour_eau)+
   tm_lines(col = "blue")+
tm_shape(detail_hydro)+
   tm_dots(col = "red")+
tm_shape(surf_hydro)+
   tm_polygons("steelblue")+
tm_shape(borders)+
   tm_borders(lwd = 2)+
   tm_basemap("OpenStreetMap")

What about history ?

The “Etat-major” map is a general map of France made, in its first version, in the 19th century. Here how to get it :

etat_major <-  get_wms_raster(shape = borders,
                            apikey = "cartes",
                            resolution = 5,
                            layer_name = "GEOGRAPHICALGRIDSYSTEMS.ETATMAJOR40",
                            filename = "etat_major")
#> 1/2 downloaded
#> 2/2 downloaded
file.remove("etat_major.tif") # I don't want to keep etat_major on disk because this is vignette
#> [1] TRUE

tm_shape(etat_major)+
   tm_rgb()+
tm_shape(borders)+
   tm_borders(lwd = 3, col = "red")+
   tm_basemap("OpenStreetMap")
#> stars object downsampled to 980 by 1020 cells. See tm_shape manual (argument raster.downsample)